Ultrasolitons: multistability and subcritical power threshold 
from higher-order Kerr terms 
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PACS 42.65.Tg - Optical solitons; nonlinear guided waves. 

PACS 42.65. Jx - Beam trapping, self-focusing and defocusing; self-phase modulation. 

PACS 42.65.Pc - Optical bistability, multistability, and switching, including local field effects. 

Abstract - We show that an optical system involving competing higher-order Kerr nonlinearities 
can support the existence of ultrasolitons, namely extremely localized modes that only appear 
above a certain threshold for the central intensity. Such new solitary waves can be produced for 
powers below the usual collapse threshold, but they can also coexist with ordinary, lower-intensity 
solitons. We derive analytical conditions for the occurrence of multistability and analyze the 
dynamics of the different kinds of fundamental eigenmodes that can be excited in these nonlinear 
systems. We also discuss the possible transitions between solitary waves belonging to different 
nonlinear regimes through the mechanism of soliton switching. 



Introduction. — A recent measurement of the instan- 
taneous higher-order Kerr (HOKE) coefficients in gases [1] 
has led to a revolutionary description of the filamentation 
of ultrashort laser pulses [2]- The filament stabilization 
is attributed to the competition of the HOKE focusing 
and defocusing contributions to the refractive index alone, 
rather than to their interplay with the ionization-induced 
plasma defocusing, which had a key role in the traditional 
interpretation [3]. This new paradigm has stimulated an 
increasing amount of work jlHll]. aimed either at testing 
the controversial results reported in Refs. [UIl], or at the 
theoretical exploration of its rich phenomenological impli- 
cations ptilT]. 

In this Letter, we will demonstrate that, within a well- 
defined parameters region, a system involving just local 
HOKE nonlinearities can support the existence of a new 
branch of localized solutions that will be called ultrasoli- 
tons. Such stationary states coexist with solitons of lower 
intensity, similar to those found for common optical me- 
dia [12]. Very remarkably, this implies the emergence of 
optical soliton multistability (OSM), i.e. the existence of 
two or more stationary states with the same power and 
different propagation constants and profiles, like in the 
systems reported in Refs. pE^ET] and in the recent work 



[22], that has appeared while we were preparing the re- 
vised version of the present paper. We will also derive 
an analytical condition on the HOKE coefficients for the 
emergence of OSM, and show that the ultrasolitons can be 
subcritical, i.e. they may exhibit powers even below the 
ordinary collapse threshold [23]. Finally, we will discuss 
the transitions among multistable states through efficient 
soliton switching processes. 

Mathematical model. — Let us consider a wave sys- 
tem evolving along the rj direction in the space of trans- 
verse coordinates ^ and Xy ^-nd assume that the complex 
wavefunction ^'(CXi^) satisfies the (dimensionless) non- 
linear Schrodinger equation 
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Fd^-n* = 0, 



where = d'^/dS,^ + d'^/dx^ and 



(1) 



(2) 
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This formalism applies to the paraxial propagation of the 
linearly polarized electric field E{x,y,z) of a laser pulse 
of mean wavenumber in vacuum fco = 27r/Ao, being Aq the 
central wavelength, in a nonlinear optical medium whose 
refractive index depends upon the intensity / = cqcE"^ 
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as n = 7io + A?! = riQ + X]g=i n2ql'' ■ In fact, motivated 
by the results of Ref. [T] for the optical response of com- 
mon gases, we will assume that the refractive index in- 
volves 4 terms of increasing powers in the intensity / of 
the beam, with alternating sign coefficients 712 , rig > and 
n4, ng < 0, that contribute to focusing and dcfocusing re- 
spectively. Wc will choose the dimensionlcss quantities 
such that /2 = A = 1- The relations between the dimen- 
sional and dimensionlcss physically-relevant quantities are 
then (^,x) = {m/\ni\Y''^{kon2){x,y), -q = (fcon|/|n4|)z, 
= {eoc\n4^\/n2Y'^E{x,y,z), An = {nl/\n4)F, 
n-G = (n|/n2)/6 and rig = (n|/n2)/8- The only free pa- 
rameters included in Eq. [l]are then /g and /g, that will 
be assumed to be positive. 

For the sake of clarity, and motivated by Refs. [Tll^fT^. 
in this Letter we will neglect the effects of multiphoton 
absorption, ionization and temporal dispersion of the op- 
tical pulses. However, we emphasize that our results may 
also be applied to multilevel atomic media, where the 
dependence of An given by Eq. [2] may be achieved by 
the coherent control of the atomic ensemble via quantum- 
engineering techniques |24| . 

Analytical condition for multistability. Let us 

now derive an analytical condition for the possible emer- 
gence of OSM. Assuming radial symmetry, we search for 
soliton solutions of the form x, v) — ^{r)e~'^'^^ , where 
r = ■\/<^^ + and jj, is the propagation constant. As 
discussed in Refs. [T^dS], /i can be identified with the 
chemical potential of an equivalent thcrmodynamical 2D 
system of = / pd^dx = J\^\^d^dx particles (for the 
optical system, N = nofcQn2p, where V is the total power 
of the optical field). Eq. [l]can then be derived by mini- 
mizing the Landau's grand potential [26] $7 = — J pd£,dx, 
where the pressure field p is 

$ I ^ 

P=-^|Vi$|2 + ^|ci>|2+^ F{U)dU. (3) 

In particular, assuming the dependence given in Eq. 
[21 the integral term in Eq. [3] can be written as 

In the case of a high-power solution with a flat-top pro- 
file of radius R , calling A — <f>(0) the amplitude of the 
large and homogeneous central region, Eq. [1] implies 

fi^-F{A'), (4) 

where we choose the arbitrary phase of the solution such 
that A is positive real. On the other hand, as shown in 
Refs. [12l[25], any such flat-top solitons obey the Young- 
Laplace (YL) equation [26], Pc — 2<t/R, where R is the 
radius of the droplet, Pc = p{0) is the central pres- 
sure and the effective surface tension can be computed 
as cr = —R~^ rp{r)dr. In the large R limit, Pc = 0, 
the gradient term in Eq. [3] can be neglected close to the 
origin, and using Eq. [4] we get 




i^(f/)df/ - A^F(A^) = 0, (5) 
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Fig. 1: (Colour on-line) Multistability domain given by Eq. 
[71 Upper inset: nonlinear refractive index correction F as a 
function of the peak intensity of the light beam Lower 
inset: detail of the region displayed within the red box depicted 
in the upper inset. 



being ^00 and Aoo the asymptotic values corresponding to 
the i? — >■ 00 droplet. 




Fig. 2: (Colour on-line) Optical power A*" vs. propagation con- 
stant for the localized nodeless stationary states of the sys- 
tem described by Eq. [H Both axes display logarithmic scales. 
The red dotted (black solid) line stands for stable ordinary 
solitons (ultrasolitons) , while the black dashed line stands for 
unstable ultrasolitons. Labeled points on the main curves refer 
to the eigenstates displayed within the insets. Insets: radial 
profiles of several eigenstates of the system, whose propaga- 
tion constants are jia = —1.46 =, jib = —0.18, /Xc = —0.20, 
Hd = —5.26, = —0.24 and /i/ = —6.67, respectively. Notice 
that all solitons depicted within each inset feature the same 
power A'^, being this a proof of the existence of OSM. 
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In particular, for media whose nonlinear response is de- 
scribed by Eq. ^ including terms up to fs, we get the 
condition 

l^lu + f^lu'^fs-^U'^O, (6) 

where U = A^. This cubic equation, giving fe and fs, 
can have either one or three real roots. The latter case will 
eventually correspond to OSM. After a long but straight- 
forward algebra, we get the following necessary and suffi- 
cient condition on /g and /§ for Eq. |6] to have three real 
solutions: 

18225/1 - 5400/2 - 77760 fefs + 20480/8 + 93312/| < 0. 

(7) 

In Fig. [U we plot the region of the {fe,fs) plane that 
satisfies such a condition, assuming fe, fs > 0. In partic- 
ular, for fe > 0.38 or fs > 0.05, condition [7] is not fulfilled 
and only one solution can be found. Conversely, for pairs 
{fe^fs) lying within the shaded region in Fig. [U three 
different branches of solitons can be found in principle, 
whose amplitude A^o for i? — >■ cx) can be obtained from 
Eq. ini In the case of oxygen (air), that was examined in 
Ref. [12], fe = 2.8 and fs = 3.9 (/g = 11.2 and fs = 34.1). 
These values fall out of the shaded region in Fig. [U and 
in fact only one branch of solitons was found in Ref. [12] . 

Hereafter, we will fix the values fe = 0.3 and /§ = 0.02, 
that lie in the OSM domain, although we have verified that 
similar results can be obtained for different choices satis- 
fying Eq. [7] The nonlinear refractive index dependence 
on the intensity of the input beam acquires then a double- 
hump structure (see insets in Fig. [1]). In fact, the emer- 
gence of OSM can be heuristically related to the appear- 
ance of the two maxima of the refractive index, suggesting 
the possible existence of families of fundamental solitons 
with limiting intensities close to the values corresponding 
to the two local maxima of F. However, the double-hump 
structure alone does not guarantee the emergence of OSM, 
as we have checked by finding values of (/g, fs) not satis- 
fying the multistability condition [7] that nevertheless lead 
to a similar double-peak structure for F. 




Fig. 3: (Colour on-line) 3D pseudocolour plots of the amplitude 
of three flat-topped beams modeled by Eq. [9] with amplitudes 
A = A'^ (left), A^ (middle), A]^ (right), after a propagation 
distance of ry = 2000, 100 and 2000, respectively. They have 
been initially perturbed with a 5% random noise. We see that 
both beams featuring and A^ remain stable while the 
one having A^^ undergoes filamentation. The spatial scales 
spanned are (5,x) £ [-100, 100]. 

With the choice fe = 0.3, fs = 0.02, we calculate the 
three roots of Eq. HI A°^^ 1.087, A^^ = 1.601, A^ = 



3.212, together with their corresponding values for the 
propagation constant, as given by Eq. S] /ij^ = —0.241, 
fi^ = —0.182, = —6.72. The superscripts o, g, u re- 
fer to the names that will be given to the three branches 
of solitons corresponding to these limiting values, namely 
ordinary solitons, ghost solitons and ultrasolitons, respec- 
tively. 

Localized stationary solutions. Ultrasolitons. 

We have solved Eq. [T] numerically, and found the localized 
stationary states shown in Fig. [5] We obtain two differ- 
ent branches of solutions that cannot be connected with 
each other, i.e., there are no bifurcations in the eigenstates' 
structure. The red dotted line stands for the ordinary soli- 
ton branch, i.e., the branch of solutions similar to those 
reported in [12] , whose lower (upper) limit corresponds to 
the critical power [^^ for self- focusing (A'^ plane wave). 
The presence of HOKE nonlinearitics gives rise to a new 
family of solutions, represented by black lines in Fig. [2l 
whose upper limit corresponds to the A^ plane wave. 
Their lower limit does not correspond to the Kerr limit, 
which means that their existence cannot be explained by 
a balance between diffraction and the leading Kerr non- 
linearity driven by /2, but rather as an interplay between 
competing HOKE nonlinearitics. To our best knowledge, 
such solutions do not have counterparts in any other non- 
linear optical system ruled by local intensity-dependent 
nonlinearities because they exist over a certain intensity 
threshold and feature both amplitudes and propagation 
constants higher (in absolute value) than those of the or- 
dinary ensemble. In other words, these solitons belong to 
a completely different nonlinear regime as compared with 
that of the ordinary branch. For all the previous reasons, 
we have called them ultrasolitons. On the other hand, 
according to the Vakhitov-Kolokolov (VK) criterium [27] 
supported by systematic simulations of propagation, we 
find that the eigenstates represented by the solid (dashed) 
line arc dynamically stable (unstable). 

The eigenstates b), c) and d) , that are represented in 
the lower-left inset of Fig. [2] feature the same optical 
power but different propagation constants and radial pro- 
files. Analogously, the flat-top eigenstates of insets e) and 
f) also correspond to the same power, thus demonstrating 
OSM even in the high power regime. Both cases demon- 
strate the emergence of OSM, generalizing to our multi- 
branch situation the definition given in [18] for a single 
continuous branch. 

On the other hand, the solution displayed in inset a) of 
Fig. [2] is the first stable ultrasoliton, and corresponds to a 
local minimum of the ultrasoliton power curve depicted in 
Fig. [2] being this a trace of the uniqueness of this solution. 
Very remarkably, this soliton features a subcritical power 
(N ~ 3.83 in dimensionless units), below the ordinary col- 
lapse threshold {N = 5.85). To our best knowledge, this 
is the first example of a stable subcritical soliton in an 
optical system with an instantaneous nonlinear response. 
Moreover, this minimum power state turns out to have 
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Fig. 4: (Colour on-line) Evolution of the peak amplitude 
A = $(0,0) of the eigenstates b) (blue dashed line), c) (green 
dashed-dotted line) and d) (red solid line) displayed in Fig. 
[21 We observe that both b) and d) fields are stable against 
small perturbations, while the eigenmode c) becomes unstable, 
decaying to the lower branch through a soliton switching mech- 
anism. Alternatively, the unstable mode c) can excite a stable 
ultrasoliton (black dotted line) by adding an initial wavefront 
curvature to the optical field. The upper (lower) row of inner 
snapshots show several pseudocolour amplitude plots of the 
optical field during the soliton switching procedure, where an 
ultrasoliton (ordinary soliton) is excited from the unstable field 
c). The square window displayed in the snapshots has a width 



the lowest possible central intensity (/ = 9.98) among the 
stable ultrasolitons, although we have found unstable ul- 
trasolitons (corresponding to the dashed curve in Fig. [2]) 
starting from central intensities as small as / = 4.88. 

In Fig. [2]we have only plotted two branches of solutions, 
the ordinary and the ultra- solitons. In fact, we have not 
been able to find numerically the third family of solitary 
waves linked to the plane wave with amplitude that 
we have theoretically predicted above. We can understand 
its non-existence by using three different arguments. 

First, we recall that, according to Refs. [121111], for large 
but finite fiat-topped solutions the central pressure Pc does 
not vanish, being compensated by the surface tension a 
like in an ordinary liquid. The actual values of a for the 
three branches can be computed as 
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0.0484 



(8) 

0.1 20i and 



Wc find: (7° = 0.0948, 
a" = 7.21. The fact that for the ghost family the surface 
tension is not real reflects the impossibility of equi- 
librating the inner pressure, which is a real magnitude, 
thus forbidding the existence of finite ghost solitons. In 
other words, the ghost solitons would not fulfill the YL 
equilibrium condition. 



Second, we have applied to the plane wave solutions 
of our system a linear stability analysis similar to that 
of Refs. [28ll29]. After a straightforward study, we have 
found that the solutions with A'^, A"^ are linearly stable, 
i.e., they do not undergo modulational instability under 
small perturbations, while the A^ lies within an instabil- 
ity window. 

Third, wc have studied numerically the propagation of 
three fiat-top beams belonging to the high power regime 
of each of the three branches. The initial condition of our 
simulations is modeled by the following function 

</> = A(0.25{[l + tanh(r + w)][l - tanh(r - w)]}), (9) 

where A = j4^, A^, A^ is the amplitude of the beam 
envelope and a; = 50 is the mean radius. We perturb the 
initial beam profiles with a 5% random noise in order to 
stimulate the onset of instability. The final states arising 
from the propagation of such initial conditions are dis- 
played in Fig. m In the left (right) picture, we show the 
surface amplitude plot of the fiat-top beam having an ini- 
tial amplitude A = A'^ {A'^), at a propagation distance 
ry = 2000. We see that such beam is stable, as it quickly 
couples to a stable eigenstate of the low-intensity (high- 
intensity) branch of ordinary solitons (ultrasolitons) dis- 
cussed above. In the middle picture, we show the outcome 
of the propagation of an initial beam having A = A^, at 
a propagation distance rj = 100. In this case, we observe 
how the spatial beam profile has been destabilized by the 
growth of the perturbations, yielding to multiple filamen- 
tation like in the Cubic-Quintic model |29| . The arising 
filaments correspond to perturbed quasi-stationary ultra- 
solitons. In this context, the onset of filamentation can be 
considered as an additional trace of the non-existence of 
the flat-top soliton with A « A^, although it would not 
be conclusive if taken alone. 

However, our present analysis is not sufficient to exclude 
that triple-stability might be found for a different choice 
of fe and /g within the multistability region of Fig. 1. 

Soliton switching. — The existence of OSM, as 
shown in Fig. 121 suggests the possibility of observing tran- 
sitions between different multistable states. Such beam- 
reshaping mechanisms may lead to soliton switching pro- 
cesses with a great potential for all-optical communica- 
tions [20]. In order to study the soliton switching in our 
system, we have simulated the free propagation of three 
perturbed solitary waves featuring the radial profiles b), 
c) and d) of Fig. 121 The optical power of these beams is 

w 30. The results of the numerical computations are 
summarized in Fig. 31 plotting the peak amplitude evo- 
lution of the eigenstates b) (blue dashed line), c) (green 
dashed-dotted line) and d) (red solid line). Even though 
wc have added a 5% random noise to their initial profiles, 
both fields b) and d) remain stable, in agreement with the 
prediction of the VK criterium applied to our system. On 
the other hand, the unstable field c) rapidly decays to a 
nonlinear mode similar to b) , even in the absence of exter- 
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nal random noise. The conversion efficiency between both 
modes, measured as the ratio between the optical powers 
of the fields, is above 90%. 

We have observed that the unstable fields always de- 
cay to the lower branch (down-switching), while the up- 
switching to the ultrasolitons realm does not occur spon- 
taneously. However, we can force such a transition by 
including a focusing quadratic phase term e"'" '^-'^'' , sim- 
ilar to that introduced by a thin lens. In this case, the 
unstable beam d) can be promoted to the upper branch 
(black dotted line in Fig. |4]) with a maximum conver- 
sion efficiency of around 30%, while the energy excess is 
radiated as a low-intensity reservoir resembling that gen- 
erated during the excitation of Townes-like waves in air 
[30] . The difference between the efficiencies of both pro- 
cesses can be understood by mode-coupling arguments. 
The modal energy transfer becomes more efficient when- 
ever both profiles and phases of the corresponding modes 
match, and the modes b) and c) of Fig. [2] have closer pro- 
files and propagation constants than the modes c) and d). 
Few stages of the soliton switching process arc highlighted 
in the insets of Fig. HI 

Conclusions. — We have shown that an optical sys- 
tem involving competing higher-order Kerr nonlinearities 
can support the existence of power multistability in the 
absence of any external potential, yielding a new class of 
solitary waves, called ultrasolitons, that may exhibit pow- 
ers below the ordinary collapse threshold. Wc have also 
proposed a mechanism of soliton switching for inducing 
transitions between different multistable nonlinear waves 
of the system, that could have potential applications in 
ultrafast optical circuits intended for all-optical communi- 
cations [20]. We hope that these results will contribute to 
stimulate the quest for a clarification of the mechanism of 
ultrashort pulses filamentation. 
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